Impurity in a Bose-Einstein condensate in a double well 



F. Mulansky, 1 J. Mumford, 1 and D. H. J. O'Dell 1 

1 Department of Physics and Astronomy, McMaster University, 
1280 Mam St. W., Hamilton, ON, L8S 4M1, Canada 

We compare and contrast the mean-field and many-body properties of a Bose-Einstein condensate 
trapped in a double well potential with a single impurity atom. The mean-field solutions display a 
rich structure of bifurcations as parameters such as the boson-impurity interaction strength and the 
tilt between the two wells are varied. In particular, we study a pitchfork bifurcation in the lowest 
mean-field stationary solution which occurs when the boson-impurity interaction exceeds a critical 
magnitude. This bifurcation, which is present for both repulsive and attractive boson-impurity 
interactions, corresponds to the spontaneous formation of an imbalance in the number of particles 
between the two wells. If the boson-impurity interaction is large, the bifurcation is associated with 
the onset of a Schrodinger cat state in the many-body ground state. We calculate the coherence and 
number fluctuations between the two wells, and also the entanglement entropy between the bosons 
and the impurity. We find that the coherence can be greatly enhanced at the bifurcation. 
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I. INTRODUCTION 

In recent years a large number of experiments have 
been performed upon atomic Bose-Einstein condensates 
(BECs) trapped in double well potentials, see e.g., p3- 
114] . Each atom can be in a superposition of being in 
both wells, allowing for studies of macroscopic quantum 
coherence [THSJ [5j IT5HT8] analogous to that found in a 
Josephson junction [H HI [13j IT9H22] . However, in con- 
trast to a traditional solid state Josephson junction, the 
microscopic hamiltonian of an atomic gas is highly con- 
trollable and is very well understood. The double well 
system is also one of the simplest ways to go beyond the 
mean-field Gross-Pitaevskii paradigm of gaseous BECs 
because the inter-well tunnelling introduces a very low 
energy scale which can easily be surpassed by the in- 
teraction energy. This tends to emphasize the discrete 
(second-quantized) aspects of the quantum state, which 
are conjugate to the phase properties, and can be seen 
in effects such as number squeezing [6, 9J. Double well 
systems are, of course, also naturally disposed to being 
used for matter- wave interferometery [21 [3J [51 [H] , with 
the eventual aim of making precision measurements. 

In this paper we consider N identical bosons trapped 
in a double well potential, and model this system via the 
two-site Bose-Hubbard model. To this we add a single 
impurity atom. The impurity can be an atom of a differ- 
ent species or internal state, but we assume that it is also 
trapped in the double well potential (but need not have 
the same tunneling rate between the wells as the bosons). 
Reference [23j considers the related, but different, case of 
an impurity atom trapped in a double well potential im- 
mersed in a uniform BEC as an example of the spin-boson 
model. In our model the bosons are limited to just two 
states (two-mode approximation), whereas in a uniform 
BEC there is a continuum of bosonic states. Another 
study which is related to ours considers an atomic quan- 
tum dot acting as a coherent single atom shuttle between 
two BECs PI]. 



One motivation for studying an impurity interacting 
with a BEC in a double well is that it can be used as 
a simple, yet concrete, model for studying the measure- 
ment problem in quantum mechanics |25j : the impurity 
is a two-state quantum system interacting with a macro- 
scopic measurement device represented by the bosons. 
The measurement device can be tuned between being 
quantum (iV small) or classical (N large). Furthermore, 
the use of a Feshbach resonance would allow the interac- 
tion between the quantum system and the measurement 
device to be tuned between being weak and strong. In 
this context we note that the interesting question of the 
classical (mean-field) limit of bosons in double well po- 
tential has been addressed in a number of papers |26l - 
|2"5| with the general conclusion that this occurs when 
N — > co, although with the caveat that quantum ef- 
fects remain important close to the separatrix j27j [28] , 
which is the boundary in phase space outside of which 
self-trapping occurs. 

A theoretical study of the same system as considered 
here has been performed by Rinck and Bruder [29] . Their 
paper, which focuses on the case of small atom num- 
ber, predicts that when the boson-impurity interaction is 
strong a tunnelling resonance occurs that involves the si- 
multaneous tunnelling of many bosons together with the 
impurity. In the presence of a tilt asymmetry between 
the two wells this resonance corresponds to the expul- 
sion of the impurity into the higher lying well. In this 
paper we extend their study to larger particle numbers 
and consider both the full many-body theory and the 
Gross-Pitaevskii mean-field theory (which is expected to 
be valid in the large N regime). We find that the Rinck- 
Brudcr tunnelling resonance is associated with a bifurca- 
tion in the mean-field theory solutions. 

Although it might seem rather fanciful to study the 
problem of a single impurity in a BEC in a double well, 
we note that a recent experiment |30j has realized an op- 
tical lattice with many bosons and one impurity per site, 
which is quite close to the situation we are considering 
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here. 



II. THE MANY-BODY HAMILTONIAN 

We begin by writing down the many-body hamiltonian 
for N identical bosons interacting with a single impurity. 
Both the bosons and the impurity are trapped in a dou- 
ble well potential (for simplicity we shall assume in later 
sections that the potential is the same for both). The 
impurity is a particle distinguishable from the bosons — 
it may be either a boson or a fermion, but its statistics 
do not matter. It can even be the same species of atom 
as the bosons, but in a different hyperfine state provided 
there is no interconversion between the hyperfine states. 
We then make the two-mode approximation (i.e. the sin- 
gle band approximation for the two-site Bose-Hubbard 
model) for the bosons and likewise for the impurity, so 
that there are four modes and 2 x (N + 1) many-body 
states in total. 

The total hamiltonian for the BEC plus impurity sys- 
tem is 



H = H A + H B + H 



BB 



H 



AB 



(1) 



where Ha is the single-particle hamiltonian for the impu- 
rity, Hb is the single-particle hamiltonian for the bosons, 
Hbb is the boson-boson interaction hamiltonian, and 
Hab is boson-impurity interaction hamiltonian. In terms 
of the field operators <&(x) for the bosons and ^>(x) for 
the impurity, we have 



H b 
Hbb 



d 3 x&(x) 
d 3 x&(x) 



2mA 

h 2 



v 2 
v 2 



V A (x) 
Vb(x) 



9b B 



2m B 

d s x&(x)&(x)$(x)®(x) 



(boson-boson interaction) 
H A b=9ab [ d 3 x&(x)&{x)${x)4>{x) 



9(x) (2) 

(3) 
(4) 

(5) 



(boson-impurity interaction) 



where gsB = ^irh^asB /mB characterizes the boson- 
boson interactions in terms of the inter-boson s-wave 
scattering length cibb and the boson mass mj, and 
cjab = 2irh 2 a,AB /m r characterizes the boson-impurity 
interactions in terms of the boson-impurity s-wave 
scattering length ciab and the reduced mass m r = 
mATHB I '(rriA with being the mass of the im- 

purity atom. The field operators in the above expres- 
sions obey the commutation relations [<&(x),ffi (x')] = 
S(x - x 1 ), [¥(x), #t( x ')] = s(x - x'), [$(x), &(x')} = 0, 
and [$(»), *(a/)] = 0. 

In this paper we shall work in the two-mode approx- 
imation. For each species the ground mode is sym- 
metric and the excited mode is antisymmetric: the en- 



ergy splitting between them is determined by tunnel- 
ing and is therefore exponentially smaller than the en- 
ergy separation to higher modes, allowing us to isolate 
just the lowest two modes. Even and odd combina- 
tions of the symmetric and antisymmetric modes give 
two new orthogonal modes which are localized, respec- 
tively, on the left and right hand sides of the double well: 
4>i{x) for the bosons and ipi(x) for the impurity, with 
/ 4>*4>jd 3 x = J ip*ipjd 3 x = Sij, where i and j can be ei- 
ther L or R. Expanding the field operators in terms of 
the localized modes we have 

*(x) = a L ip L (x) + a R ip R {x) (6) 

$(x) = b L <pL{x) + b R <t>R{x), (7) 

where the mode operators for the bosons obey the com- 
mutation relations = 5^, and [Sj,5j] = 0. A simi- 
lar set of commutation relations hold for the impurity op- 
erators, and the boson and impurity operators commute 
with each other. Substituting these expansions into the 
hamiltonian we obtain the two-mode hamiltonian, which 
can be written, up to constant terms, as [29] 

U -2 
H = -AN - JB - J a A 
4 

W - Ae - Ae a 

+ — ANAM+—AN+—AM (8) 

where we have defined: AM = a R a R — o^cil as the num- 
ber difference operator between the two wells for the im- 
purity (its eigenvalues are ±1); AN = b R b R — b\b L ditto 
for the bosons (its eigenvalues range from — N to N in 
steps of two); and A = a L a R +a R a,L and B = b L b R +b R i)L 
as the hopping operators for the impurity and bosons, re- 
spectively. The parameters that appear in the two-mode 
hamiltonian are defined as 



Ul,r = 9bb 



d 3 x\<j> L , R \ i 



W L . 



B 



9ab / d 3 x\(p LtR \ 2 \^ LiR \ 2 



J = - J d'xcpl x 
Ae = e R - e L . 



2m i 



Vb(x) 



(9) 
(10) 

(12) 

Ul/r is the intra-well interaction energy for the bosons, 
and Wl/ r is the intra-well interaction energy between 
the bosons and the impurity. For the rest of this paper we 
shall assume that Ul — U R = U and Wl = W R = W. J 
and Ae are, respectively, the bosonic hopping energy and 
difference in zero-point single-particle energies between 
the two wells. J a and Ae a are the equivalent quantities 
for the impurity. The single-particle energy c.l,r for the 
bosons is defined as 



£l,r 



d 3 x(f> 



L,R 



2m, 



-V 2 + V B (x) 



->L,R ■ 



(13) 



The quantity Ae can be regarded as an imbalance or tilt 
between the two wells arising, e.g., due to gravity if one 
well is lower than the other. 
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Note that the expression ^ for the two-mode hamilto- 
nian neglects small cross terms such as those that depend 
on integrals like J d?xcb* L 4>L , 4 ! L' i l } R- Furthermore, we shall 
also assume that the trapping provided by the wells is 
tight enough that the single-particle energies c-l,r dom- 
inate the interaction energies. This allows us to neglect 
changes in the mode wave functions, and hence changes 
in the parameters ([£])-( 12 ), as the particle number in each 
well varies. 

Let us find the matrix elements of H in the Fock basis 
\Mr,Nr), where Mr is the eigenvalue of the impurity 
number operator for the right well Mr = a R a,R, and 
Nr is the equivalent eigenvalue for the bosons. Due to 
number conservation, we only need to specify the number 
of particles in the right hand well. We find that the 
general matrix element is 



(M' R N' R \ 



-Ac 
,2 



AN 



H \M R N R 
U AN 2 



-AM 



W 



AN AM) 6 



N'NrOM'Mr 



Ac 

4~" ~2 ~"~ ' 2 

-Jy/N R (N -N R + 1)6 n ' r n r -i5m' r m r 
-J^(N-Nr)(Nr + 1)6 n 

' R N R + l&M' R M R 

—J a \j M R [M -M R + 1)8 n > rNr 8m> h Mr-i 
-J a ^/(M - Mr)(Mr + 1)5 n , Nr 5 m , Mr+1 



(14) 



Ne 



where M = 1 is the number of impurities, AN 
N L and AM = Mr - M L . 

Diagonalizing the matrix specified by Eq. ( 14 ) gives 



the eigenvalues and eigenvectors of the many-body two- 
mode hamiltonian. In Figure[l]we plot the eigenvalues as 
a function of the tilt Ae for the case of six bosons for both 
positive and negative values of the boson-impurity inter- 
action W. The avoided crossings correspond to places 
where a particle hops from one well to another as the tilt 
is changed. Similar pictures have previously been made 
for the BEC-impurity system by Rinck and Bruder [25] . 
who noted that when the boson-impurity interaction W 
is repulsive enough the impurity can be forced to tunnel 
uphill against the gradient set by Ae a , i.e. the impurity 
can be expelled from the BEC into the higher lying well. 

The problem of a BEC in a double well potential can 
be mapped onto a pendulum model, or, equivalently, a 
single particle in a periodic potential. As explained in 
[2"T] . in this latter picture the tilt Ae plays the role of 
the quasimomentum of the particle, and the eigenvalue 
structure seen in Fig. [I] as a function of tilt can then 
be viewed as a band structure plotted as a function of 
quasimomentum. The energy level structure shown in 
Fig. [I] does not have the usual periodic form we expect of 
a band structure, but it can be made periodic by applying 
a simple transformation |27j . 
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Figure 1. (Color online) Eigenvalues of the many-body hamil- 
tonian as a function of the tilt Ae (= Ae a ) for TV = 6 bosons 
and a single impurity. Each panel has a different value of the 
boson-impurity coupling: (a) W/U = —4, (b) W/U = —2, (c) 
W/U = 0, (d) W/U = 2, (e) W/U = 4. The values of the 
other parameters are J/U = J a /U = 1.5. 



III. MEAN-FIELD APPROXIMATION 



Let us now use the mean-field approximation to calcu- 
late the equations of motion for the number and phase 
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differences of the condensate and impurity wave functions 
between the two wells. The stationary solutions of these 
equations provide natural quantities to compare with the 
eigenvalue structure of the second-quantized theory be- 
cause both quantities are stationary during time evolu- 
tion. The relationship between the second quantized the- 
ory presented in Section [TTJ and the mean-field theory 
presented in this section, is analogous to that between 
quantum and classical mechanics for a single particle: 
the N — > oo limit of the many-body theory gives the 
mean-field theory, and is equivalent to the H — > limit of 
single particle quantum theory [55]. In fact, the precise 
relationship is NJ/U = {S/H) 2 where S is the classical 
action [27] . 

In the presence of Bose-Einstein condensation the 
bosonic field operator can be written = ^o^) + 

5Q(x), where &o(x) is the so-called condensate wave func- 
tion. A description of the bosons solely in terms of <E>o 
corresponds to a mean-field approximation. Making this 
replacement in the hamiltonian gives the mean-field 
energy functional H (because there is only one impu- 
rity atom, we can also replace its field operator by its 
wave function '5 — > \P without involving any approxi- 
mation). The equations of motion for the impurity and 
condensate wave functions can be found by taking func- 
tional derivatives of H, i.e. ihd^/dt = SH/5^* and 
ihd^o/dt = (5iJ/(5$Q- The result is a Schrodinger equa- 
tion for the impurity 



„0¥ 
lH -dt = 



2mA 



V A (x)+g A B |* r 



(15) 



coupled to a Gross-Pitaevskii equation for the bosons 



ih 



dt 



2m B 



V 2 + V B (x)+g B \%\ 2 + g AB m 



(16) 

In order to make the two-mode approximation in the 
mean-field case, we replace the mode operators a L / R and 
b L / R that appear in Eqns. (JgJ) and ^ by the complex 
numbers cll/r an d &l/_r 



a L/R 



= JM. 



?L/R 



N. 



l/rc 



(17) 
(18) 



Substituting these forms into the Schrodinger and Gross- 
Pitaevskii equations we obtain the equations of motion 



Ae a W AJ a Fcosa 
a = + 2 ~T Z + -£ / 

h h h y/l - 4Y 2 

Y = — \/l - AY 2 sina 

h 

■ Ae U „ „W„ AJ Z cos/3 
P= it +2-Z+2— Y+ - — H 21 

h h h h y/N 2 - 4Z 2 



(19) 
(20) 



Z = --y/N 2 ~AZ 2 sin/3 

h 



(22) 



where we have defined the variables a = — olr, Y = 
AM/2, (3 = f3 L - (3 R , and Z = AN/2. Recognizing that 



the canonically conjugate pairs of variables are {a, HY} 
and {(3,hZ}, the equations of motion can be expressed 
in the form of Hamilton's equations 



1 OH 

hW 

1 OH 

YdZ 



Y = - 



Z =- 



1 dH 

h da 
1 dH 

h~dfi 



(23) 
(24) 



where H in terms of the new variables is 

H = UZ 2 - J\/N 2 - AZ 2 cos/3 - J a \/l - AY 2 cos a 
+2WYZ + Ae Z + Ae a Y . (25) 

This mean-field hamiltonian function can be compared 
term by term with the second-quantized version given in 
Eq. @. 

In order to gain some intuition, let us first consider the 
case where W — 0, so that the impurity and the BEC are 
not coupled to each other. The impurity is an elementary 
two state system, analogous to, e.g., a spin or a two-level 
atom. The equations of motion for the impurity can be 
solved exactly to give 

Y(t) = Y sin (wi 

) (26) 
-2Y cos(w imp i + <f>a) 



a(t) = arcsin 



yj 1 - AY 2 sin 2 (w imp i + 4> a ) 



(27) 



where the constant (j) a gives the initial phase of the 
motion, the constant Yq sets the magnitude of the os- 
cillations of the amplitude Y(t), and must lie in the 
range —1/2 < Y < 1/2, and the angular frequency 
^imp = 2J a /h of the oscillation is given by the bare hop- 
ping frequency. For simplicity we have assumed that the 
tilt Ae a /h = 0. 

When W ^ 0, the coupled motion of the impurity and 
the BEC can be complicated (indeed, it is expected to 
be chaotic — see later). A simplified situation arises if 
the BEC is static, i.e. if J = so that its tunneling is 
switched off and Z is locked at a particular value. We 
may then ask, how does the presence of the BEC affect 
the tunneling frequency of the impurity? The answer is, 
not very much. As can be seen by inspection of Eqns. 
( 19 )-( 22 ) , the effect of a static BEC upon the impurity 



is exactly the same as the tilt term containing Ae a , which 
does not affect the frequency of motion. In fact, if Z = 
the impurity is completely unaffected by the BEC. Thus, 
the impurity does not acquire an 'effective mass' through 
its interaction with the BEC. This result probably only 
holds in the Bose-Hubbard limit considered in this pa- 
per. When the mode wave functions are allowed to be 
modified by interactions it is likely that the impurity will 
acquire an 'effective mass' which will affect Wi mp . 

Let us now turn to the BEC component. In the absence 
of the impurity, the equations of motion (21 1 and ( 22 ) for 



the BEC correspond to the celebrated Josephson equa- 
tions |21U31] known from the theory of superconductivity. 
Providing the atom number imbalance between the two 
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wells is much smaller than the total number N (low en- 
ergy regime) , and the inequality NU >• J is obeyed [21] , 



we can replace y/N 2 
by N to leave 



AZ 2 in the mean-field hamiltonian 



H BEC « UZ 2 - J Ncos (i. 

For small angles, this yields the equations of motion 

2U „ 




Z 



JN 



(28) 

(29) 
(30) 



whose solutions arc 



Z(t) = Z Q sin (u;pi as £ + 4> z ) 
f3(t) = p sin (uj pias t + <f>p) 



(31) 
(32) 



The constants /?o, Z$ and 4>z are set by the boundary 
conditions, and the amplitude of the number difference 
oscillations must lie in the range N/2 < Zq < N/2. The 
frequency of the oscillations predicted by the above equa- 
tions of motion is 



Wplas — 



V2JUN 



(33) 



and is known as the plasma frequency. We see that, due 
to the interactions between the bosons, the frequency of 
their oscillation depends on their total number N and 
hence can, counter-intuitively, be much faster than that 
of the impurity (assuming all the parameters for the im- 
purity are chosen to be equal to their counterparts for 
the bosons). 



The hamiltonian ( 28 1 for the BEC has the same form 



as that of a pendulum 



-^pend 



2ml 2 



mgl cos P 



(34) 



where m is the mass of the bob, I the length of the (mass- 
less) rod, and g is the acceleration due to gravity. For the 
case of the pendulum, /3 becomes the angular displace- 
ment from the downward vertical, and it is conjugate 
to the angular momentum p. The analogy between the 
pendulum and the BEC is completed by relating p to the 
population imbalance via p = hZ. 

From the position of the factor N in the 'potential en- 



ergy' term in Eq. (28), it is tempting to associate it with 



the length I of the pendulum, so that I oc N [H [T!J] . How- 
ever, this is wrong because it implies that the frequency 
Wpend = ygjl of the pendulum reduces with increasing 
N, in contradiction to what w p i as predicts. Because N 
is dimcnsionlcss, we can instead let N = L/l, where L 
is a constant with units of length, e.g. the length corre- 
sponding to N=l. Then the length of the effective pen- 
dulum is inversely related to the number of atoms and 
Wpend = \J Ng/L, as required. A similar approach can 



be applied to the more general case of Eq. (25 1, where 
the terms \J N 2 — AZ 2 imply a pendulum whose length 



changes during the motion [TS] . The length of this pen- 
dulum increases with Z. 

The above considerations lead us to a picture for the 
combined BEC-impurity system as two coupled pendula 
with natural frequencies ui\ = \J g/l and u>2 = y/g/L, 
where I oc 1/JV. When N > 1 the length of the BEC 
pendulum is much shorter than that of the impurity pen- 
dulum so that I <C L. The coupling 2WYZ between the 
two pendula depends upon the product of the two an- 
gular momenta. This implies that we cannot think of 
two pendula coupled by a spring, because that would 
lead to a coupling term that depends on the difference in 
the angles a — (3. An alternative model which does have 
a coupling of the correct form is the double pendulum, 
where one pendulum is suspended from the other. The 
hamiltonian for the double pendulum is derived in the 
Appendix and can be written (when mi 7712) 



-Hdp = 



9 
Pi 



2 \ rmll 
—m\gl\ COS' 



2 

P'l 



Pi P2 



m 2 l 2 



2 n»i h h 
m 2 gl2 cos 8 2 



cos 



(35) 



where p\ and P2 are the angular momenta conjugate to 
the angles 8\ and 6*2- For small deviations from the sta- 
tionary solutions the cos(£?i — 82) term is approximately 
constant. For example, when both pendula are point- 
ing down 81 = 6*2 = 0, so that cos(#i — 8 2 ) = 1. When 
one pendulum is pointing down and the other is point- 
ing up cos(6'i — 62) = — 1. In these situations the double 
pendulum model serves as a qualitative model which can 
help guide our intuition for the BEC-impurity system. 
Indeed, a well known feature of the double pendulum is 
that it can display chaotic motion, and we expect this to 
also be true of the BEC-impurity system in the mean- 
field regime. However, one shortcoming of the double- 
pendulum model is that one cannot choose the strength 
of the coupling independently of the properties of the 
individual pendula. 



IV. STATIC SOLUTIONS TO THE 
MEAN-FIELD EQUATIONS: SWALLOWTAIL 
LOOPS AND PITCHFORK BIFURCATIONS 

In this section we compute the static solutions to the 
mean-field equations of motion. These can be found by 
setting the left hand sides of Eqns. ( 19 )— ([22]) equal to 
zero. We shall ignore the rather specialized solutions 
Z = ±N/2 and Y = ±1 to Eqns. (|20| and ((22]), and con- 
sider only the solutions a = {0, tt}, and f3 = {0, ir}. This 



leads to four different combinations for a and /? which 
we can insert into the other two equations of motion, and 
solve numerically. Finally, the energies of the static so- 
lutions are found by plugging them into the hamiltonian 
(25 ). The results of this procedure are plotted in Fig. [2] 
as a function of the tilt Ae (we have set Ae a = Ae for 
simplicity). Each panel corresponds to a different value 
of the boson-impurity interaction W, and for comparison 
we have also included the many-body eigenvalues for a 



6 




1. 

0.5 


-0.5 





(c) 







-5 5 
Ae (units of U) 

Figure 2. (Color online) Energies of the static solutions to 
the mean-field equations 122| as a function of the tilt 

Ae (= Ae a ). The various solutions are characterized by their 
phase differences: a — ft = (black circles); a = it, f3 — 
(orange squares); a = 0, /3 = it (blue diamonds); and a = 
P = 7r (red triangles). Each panel has a different value of the 
boson-impurity coupling: (a) W/U = —4, (b) = —2, 

(c) W/U = 0, (d) W/U = 2, (e) W/(7 = 4. All panels 
have J/U = J a /U — 1.5. We have also included the many- 
body eigenvalues (thin lines) for N = 6 bosons and a single 
impurity. 
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Figure 3. (Color online) A zoom-in of Fig. |5] showing the birth 
of the loop in the lowest band as the boson-impurity interac- 
tion strength W is varied: (a) W/U = 2.0; (b) W/U = 2.122; 
(c) W/U=2.2. The solid lines are given by solving F(6) = 0, 
as given in Eq. (40 1, and the dots are given by solving the 3rd 
order Taylor expansion Fs(9), as given in Eq. (41 1. The other 
parameters are set at N = 6, J = J a = 1.5 U. 




system with N = 6 bosons. Even for this small number 
of particles, the lowest and highest mean-field solutions 
cling closely to the lowest and highest many-body en- 
ergies. The two intermediate mean-field solutions (and 
also the loops in the lowest and highest mean-field so- 
lutions) do not each cling to a single many-body state, 
but rather, they pass right through avoided crossings, 
jumping between two many-body states in the process. 
This behaviour can also be seen in Figs. [6] and [7J The 
avoided crossings may be viewed as a tunneling effect: 
when either of the relative phases a or f3 is close to it 
then we are near the barrier top of a cosine potential [see 
the hamiltonian Eq. (25)] and tunneling corrections to 



the mean- field solution become important. 

The mean-field solution with a — j3 — (plotted with 
black dots) corresponds to both pendula pointing down, 
and has the lowest energy. The next solution up for our 
parameters has a = 7r and j3 = (orange squares), and 
corresponds to the BEC pendulum pointing down and 
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the impurity pendulum pointing up. The solution with 
a = and /3 — tt (blue diamonds) corresponds to the 
BEC pendulum pointing up and the impurity pendulum 
pointing down, and is the third most energetic solution 
for Ae = (various branches of different solutions can 
cross for Ae ^ 0). The highest energy solution has a — 
/3 = t: (red triangles), corresponding to both pendula 
pointing up. 

The structure of the mean-field solutions is quite rich 
and changes qualitatively as the parameters change. As 
can be seen in Fig. [2j the bands can contain swallowtail 
loops. These are a manifestation of the nonlinearity of 
the mean-field theory and have been studied previously 
in the context of BECs in double well potentials [32] , in 
the band structure of BECs in optical lattices [33M36| . 
and also in the band structure of non-interacting atoms 
in optical cavities [37]. A recent experiment on a BEC 
trapped in a double well potential has seen evidence for 
the existence of the loop structure through the violation 
of adiabaticity [38] . 

Of particular significance is the fact that the lowest 
band (black dots in Fig. [2]) undergoes a bifurcation and 
develops a swallowtail loop as the magnitude of W is 
increased beyond a critical value, ±W C (i.e. for both at- 
tractive and repulsive boson-impurity interactions). This 
suggests that the ground state of the system changes sig- 
nificantly at W c . A zoom-in of the birth of the loop at 
the birfurcation is shown in Fig. [3j In the top panel (a) 
we have W < W c , and the band is a smooth curve. In the 
middle panel (b) we have W = W c , at which point a cusp 
forms at zero tilt, heralding the birth of the loop. In the 
bottom panel (c) we have W > W c , and the band con- 
tains a loop. For the case of a pure BEC (no impurity) 
in a double well potential, a loop can also appear in the 
lowest band, but only when the interboson interactions 
are attractive (U < 0). Here it can occur for entirely 



repulsive interactions (both U > and W > 0). 

Let us now focus on the lowest band, which is defined 
by the phase differences a = j3 = 0. The lowest band is 
given by the simultaneous solution of the equations 



Ae a -I- 2WZ + 4 J a 
Ae + 2UZ + 2WY + 4 J 



Y 



VI - 4F 2 
Z 



^N 2 - 4Z 2 



= (36) 
= (37) 



To remove the square roots we introduce the variables 9 
and d>: 



Y = sin 6 
2 Y 

N 

Z = — sin ( 
2 



n 1 




n 




< 






2 


7T 




7T 




9 < 




2< 




2 



(38) 
(39) 



Eqns. (36) and (37) can then be combined to give the 
function 



F(6) =Ae + NU sin 9 + 2 J tan 6 

Ae a +NWsin9 

- W 



= = 

1 M aZ + {ke a + NW sin9) 2 
where we have made use of the ide ntity sin 
tan</>/-\/l 



(40) 



F{9) = 0, 



tan 2 <f>. Solutions of Eq. (|40j), 
are plotted as solid curves in Fig. [3] 

The lowest band typically has small population differ- 
ences Z between the two wells, and so we can restrict our 
attention to small values of 9. A plot of F(0) versus 9 
reveals that for small 9 it is either linear, or for some pa- 
rameter values it can develop a cubic structure. This fits 
in with what we see in Fig. [3] because there is either one 
or three solutions at each value of the tilt Ae, depending 
upon whether W is less or greater than W c . To obtain 
the explicit cubic equation describing the loop, we make 
a Taylor expansion of F{9) up to third order about 9 = 
to give 



W) =g I U-NU 



U a2 NW 2 ( (Ae a2 + 4J a2 ) + 12N 2 W 2 ( J a2 - Ae a 



(Ae« 2 +4J" 2 ) 7/2 




Ae a J a2 N 2 W 3 



(Ae a2 + 4J a '' 



(41) 



,5/2 



2J + NU- 



U aZ NW 2 
(Ae a2 +4J a2 



3/2 



9 + Ae 



Ae a W 



\/Ae Q2 + 4J a2 



The solutions to the equation F 3 (9) = are plotted as 
the red dots in Fig. [3j As can be seen, there is good 
agreement with the solutions of the full function F(9) = 
0. 

Let us calculate the critical value of the boson-impurity 
interaction W c at which the bifurcation in the lowest 
band occurs. From Fig. [3j we see that the loop is born 
at Ae = Ae a = 0. Furthermore, consideration of the 



way F(9) evolves between a linear and a cubic function 
shows that the bifurcation occurs when the first deriva- 
tive of F{9) or F 3 (9) vanishes at 9 = 0. Thus, we find 

W c l2J a (2J ~\~ 

W = \Jnu{u +n )- (42) 

For example, when J = J a = 1.5 U, and N = 6 we find 
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Figure 4. (Color online) The atom number difference between 
the right and left hand wells for the lowest band plotted as a 
function of the tilt Ae (= Ae a ). The left hand column gives 
relative atom number difference Z/N = (Nr — Nl)/2N for 
the BEC, and the right hand column gives the probability 
difference Y = (Mr — Ml)/ 2 for the impurity. Each row has 
a different value of the boson-impurity interaction W, and 
corresponds to the equivalent panel of Fig. [3] (a) W/U = 2.0; 
(b) W/U = 2.122; (c) W/U=2.2. The other parameters are 
set at N = 6, J = J a = 1.517. 



W c RJ 2.121 U, which is the value used in panel (b) of 
Figs. H and g) 

What value does the ratio W c /U take in a realistic ex- 
perimental situation? Let us consider the two key exper- 
iments, Gati et al 4 J and Levy et al [8], where macro- 
scopic tunneling and allied Josephson effects were first 
seen in single bosonic Josephson junctions. In the for- 
mer experiment N = 1150 and J/NU — 1/30, and in 
the latter experiment N = 10 5 and J/NU = 1/600. Al- 
though there was no impurity present in either of these 
experiments, if for the sake of argument we assume that 
J a = J, then we find that W c /U = 9, and W c /U = 18, 
respectively. If it doesn't occur naturally, this factor of 
10 between the boson-impurity and the boson-boson in- 
teraction energies can be achieved using a Feshbach res- 
onance, see, for example, the experiment [39 on the in- 
ternal Josephson effect which employed a Feshbach reso- 
nance to control interactions between different hyperfine 
states of Bose-Einstein condensed 87 Rb atoms. 

The lowest band is defined by the condition that the 
relative phases between the two wells are zero. However, 
the relative number differences Y and Z between the two 
wells can vary along the band. It is particularly inter- 
esting to see how the number differences behave in the 
vicinity of the loop in the lowest band, and this is il- 
lustrated in Figs. [4] and [5] Consider Fig. [4] first, which 
shows how Z and Y vary with the tilt Ae. The top panel 
(a) has W < W c , and there is only a single solution for 
Y and Z at each value of Ae, as expected. Because we 
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Figure 5. (Color online) The atom number difference be- 
tween the right and left hand wells for the lowest band plotted 
as a function of the boson-impurity interaction strength W . 
The left hand column gives relative atom number difference 
Z/N = (N R - N L )/2N for the BEC, and the right hand col- 
umn gives the probability difference Y = (Mr — M£)/2 for 
the impurity. The top row has zero tilt: Ae = Ae a = 0, and 
the bottom row has a finite tilt: Ae = Ae a = U. The other 
parameters are set at iV = 6, J = J a = 1.5 U. 



have chosen a positive value of W, which corresponds 
to repulsive boson-impurity interactions, the BEC and 
the impurity have opposite dependences upon Ae (we 
have set Ae a = Ae for simplicity). For our parameters, 
the BEC has the greater probability of occupying the 
lower well, so that Z = (Nr — Nl)/2 is positive when 
Ae = 6r — cl is negative. Conversely, the impurity has 
the greater probability of occupying the upper well, and 
so Y = (Mb, — Ml)/2 is negative when Ae is negative. 
The middle panel (b) has W = W c , and features a very 
sudden change in Y and Z at zero tilt. Finally, the bot- 
tom panel (c) has W > W c and the population difference 
has developed a portion which is folded back to give the 
classic "S" shape associated with a fold catastrophe [37] . 
This latter structure can generally be expected to lead 
to a hysteresis effect when the tilt is swept through zero 
in one direction versus the other direction. Consider, for 
example, the case when W > W c , and Ae is large and 
negative. Then Z will begin on its upper branch and Y 
on its lower branch. Increasing the tilt, both will follow 
their respective branches until each branch vanishes (i.e. 
each curve folds back to form the middle branches) at a 
finite positive value of the tilt. The system is then forced 
to make a jump to the lower and upper branches, respec- 
tively (or even to other bands). Conversely, if the sweep 
is performed from positive to negative value of the tilt the 
jump will occur at a finite negative value of the tilt. This 
sudden disruption to the evolution of the populations as 
the tilt is changed means that adiabatic evolution is im- 
possible in the presence of loops [331 ISZl |35J HO] , at least 
in the mean-field approximation. 

Whereas Fig.|4]shows the dependence of the population 
difference upon the tilt, Fig.[5]shows its dependence upon 
the boson-impurity interaction W. We see that the latter 
situation gives rise to a pitchfork bifurcation, as expected 
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from the cubic form of F{9). In fact, because Fig. [5] 
shows a range of W that extends from positive to negative 
values, we find back-to-back pitchfork bifurcations due to 
the fact that there are bifurcations at ±W C . The top row 
of Fig. [5] is for zero tilt, and gives a symmetric pitchfork, 
whereas the bottom row is for a finite value of the tilt 
and gives a broken pitchfork. The bottom row shows that 
there can be qualitative differences between the behavior 
of the BEC and the impurity as W is swept through 
zero (if the tilt is non-zero). For example, if the BEC 
starts off on the lowest branch for W < 0, and W is 
swept through to W > 0, it remains on the lowest branch 
(adiabatic evolution) and Z is negative throughout, i.e. 
the BEC remains in the lower-lying left well. Meanwhile, 
if the impurity also starts off on its lowest branch for 
W < (as it very well might do if the BEC starts off on 
its lowest branch, because then it is in contact with the 
BEC, thereby lowering the energy), and W is then swept 
to positive values, we find that the impurity can also 
evolve adiabatically, but is transferred into the higher 
lying right well by the end of the sweep. 

The top row of panels in Fig. [5] which is for zero tilt, 
is perhaps even more interesting because it illustrates a 
spontaneous breaking of the left/right symmetry of the 
number differences Y and Z when \W\ > W c . In the 
many-body case the system can in principle be in a su- 
perposition of more than one branch of the pitchfork, 
but in the mean-field description it must choose just one 
value for each of the number differences. Note that when 
|W| < W c we expect both Y and Z to occupy the equiv- 
alent branch of their pitchforks, because that lowers the 
energy. Conversely, when \W\ > W c they will occupy 
opposing branches, because again this lowers the energy. 
The middle branch is expected to be unstable. 

A pitchfork bifurcation of the atom number difference 
between two modes in a BEC has recently been exper- 
imentally investigated by the Oberthaler group [42]. In 
their case, the bifurcation was at the transition between 



the Rabi and Josephson regimes (see Section VI I of the 
internal Josephson effect where, rather than a physical 
double well potential, the two modes are provided by 
two spin states. 

In Fig. [6] we have plotted the many-body and mean- 
field energies as a function of the boson-impurity inter- 
action W. To sec how this figure fits in with the figures 
showing the swallowtail loops, imagine adding a third 
axis to Fig. [2] that comes out of the page and along which 
W is increased. Then Fig. [6] shows a slice through this 
3D space for a fixed value of Ae(= Ae a ). The loop in the 
lowest band (black circles) now presents itself as a pitch- 
fork bifurcation (which is broken due to the finite value 
of Ae chosen in Fig. |6| . The smaller dots are the other 
mean-field bands (the color and symbol scheme follows 
that of Fig. [2]). For comparison we have also included the 
many-body eigenvalues (solid curves) . To avoid confusion 
wc remind the reader that the hamiltonian QSj> neglects 
constant factors that only depend on the total number 
of atoms. This explains why the ground state energy in 
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Figure 6. (Color online) Dependence of the many-body and 
mean-field energies upon the boson-impurity interaction W 
for a finite value of the tilt. The solid curves are many-body 
eigenvalues, and the dotted curves are the static solutions to 
the mean-field equations: a — /3 = (black circles); a = ir, 
P = (orange squares); a = 0, /3 = n (blue diamonds). This 
figure only shows the low energies and so the a — /3 = n 
solution does not appear. The other parameters are N = 30, 
J = J a = 1.5 (7, Ae = Ae a = 717. 



Fig.|6]decreases as W increases. The true boson-impurity 
interaction energy (TV + ANAM)W/2 tends to zero for 
large W, but here we have dropped the NW/2 part. 

As can be seen in Fig. [6j the mean-field static solution 
we have been referring to as the "lowest band" (defined 
by a — f3 = 0) is not always the lowest in energy: for 
certain values of W branches of the solution with a = 
and /3 = 7r (blue dots) can have a lower energy than 
the top and middle branches of the a = p = solution. 
Nevertheless, the lower branch of the a = (3 = solution 
still has the lowest energy, and clings very closely to the 
many-body ground state energy. 

It is notable that the middle branch of the pitchfork 
in Fig. [6] forms a lower envelope bounding a region of 
avoided crossings (with the exception of the two small 
regions where the a — 0, (3 = tt band pierces the pitch- 
fork) . This is in accordance with studies [32] of BECs in 
double well potentials with no impurity, where the mean- 
field solutions bound regions of avoided crossings. 

In order to gain a detailed understanding of the behav- 
ior of the system as a function of W, we have plotted the 
N = 6 case in Fig. [7] We have reduced the value of J/U 
from that used in Fig. [6] so that the eigenstates of the 
hamiltonian correspond quite closely to the eigenstates 
of the number difference operators AN and AM except 
close to the avoided crossings (note that all apparent 
crossings are in fact avoided crossings), which makes the 
picture easier to interpret. For large values of W the ener- 
gies are dominated by the (W/2) AN AM boson-impurity 
interaction. This gives straight lines as a function of W 
whose gradients are AN AM/2, and this fact allows us to 
identify which asymptotic state is which. When N = 6 
we only have the possibilities AJV = ±6, ±4, ±2, 0, and 
AM = ±1. This is what we see in Fig. [7] where pairs of 
states can be identified at large W with gradients — 3W, 



10 



0.0- 




W (units of U) 

Figure 7. (Color online) Dependence of the many-body and 
mean-field energies upon the boson-impurity interaction W 
for a finite value of the tilt Ae = Ae a — 7U. This picture 
differs from Fig. [6] in that we have far fewer atoms (N = 6), 
and we have decreased the ratio of J to J = J a = 0.5 U. 
We have also plotted the energies of the relevant Fock states 
(dashed lines). The solid curves, which mostly lie on top of 
the Fock state energies, are the eigenvalues of the many-body 
hamiltonian. The black circles give the a — ft = static 
solutions to the mean-field equations (for simplicity we have 
not shown the other branches). 

—2W, —W, and zero (we only show the lower portion of 
the energy space and so the states with positive gradi- 
ents at large W are not shown) . Consider the two number 
states with gradient — 3W; the lower energy state of the 
pair has all six bosons in the lower (left) well so that 
AiV = —6 and the impurity in the higher well so that 
AM = 1. The higher energy state of the pair has the 
converse. For large enough W these two states become 
the ground state and the first excited state. Fig. [7] allows 
us to connect mean-field solutions to many-body states: 
the lowest branch of the pitchfork accurately gives the 
ground state energy of the hamiltonian for all W; the 
middle branch of the pitchfork gives the number state 
which has the same magnitudes of A7V and AM as the 
ground state, but with reversed signs, and which for large 
W (i.e. beyond the last avoided crossing at W as 33 U) 
becomes the first excited state of the hamiltonian, and 
the top branch gives the number state with A7V = and 
A M — — 1. When the tilt is zero the middle and lower 
branches of the pitchfork become degenerate with each 
other. 



V. MANY-BODY EIGENSTATES 

The many-body eigenstates can be written as super- 
positions of Fock states \Mr,Nr) 

\* j )= E C m r ,n r Wr,N r ), (43) 

M R .N R 

where the index j labels the j th eigenstate. For every 
possible value of Nr for the bosons, there are two possible 




Figure 8. (Color online) The first four many-body eigenstates 
for W = 0. The ground state is shown in panel (a), the first 
excited state in panel (b) and so on. The red circles/blue 
triangles are the amplitudes for the impurity to be in the 
left/right well. The total number of atoms is N = 100 and 
Nr is the number of atoms in the right well. The hopping 
energy is set to J = J a = 10 U and the tilt is zero. Accord- 
ing to hamiltonian Q, the energies of these four states are 
-0.0995, -0.0975, -0.0947, and -0.0927 in units of UN 2 . 



values of Mr, corresponding to the impurity being in 
either the left or the right well. Some examples of these 
eigenstates are shown in Figs. [H] and [9] The red circles 
correspond to the impurity being in the left well, i.e. 
Co N R -> an d the blue triangles to it being in the right well, 
i.e. C\^m r - In Fig. [8] the boson-impurity interaction W = 
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Figure 9. (Color online) The ground state for W/U = — 10 
(a), and W/f/ = 10 (b). The red circles/blue triangles are 
the amplitudes for the impurity to be in the left/right well. 
The total number of atoms is N — 100 and Nr is the number 
of atoms in the right well. The hopping energy is set to J = 
J" = 10 U and the tilt is zero. According to hamiltonian (TsJ> , 
both the positive and negative values of W lead to the same 
ground state energy of —0.101 in units of UN 2 . 



and so in the ground state shown in panel (A) the circles 
and triangles sit on top of each other. Although the state 
is peaked around AN = 0, where AA~ = Nr — Nl = 
2Nr — N, the many-body state has a finite width in AN, 
unlike the mean-field state. This width can be estimated 
by observing that the low-lying states of the pendulum 
model discussed in Section [TTl| see an essentially harmonic 
potential. Applying the results of the quantum harmonic 
oscillator gives the ground state wave function (for either 
one of the impurity states, and assuming W — 0) as 



1 



2U 



n 1 / 4 \JN 



exp 



2? 
2 



2U 

In 



(44) 



where Z = AN/2. In the case of a finite number of 
bosons this wave function is sampled discretely, but nev- 
ertheless it gives a good approximation to the exact result 
shown in Fig. [8] (A) when N is large. In particular, the 
width varies as (JN/2U) 1 / 4 . In the interaction domi- 
nated regime J <C U, the probability distribution for the 
ground state becomes narrow and tends to a single Fock 
state for which there are exactly N/2 bosons in each well. 

The symmetry of the hamiltonian Q means that its 
eigenstates have definite parity. The ground state shown 
in Fig. [8] (a) is an even function of the impurity location, 
whereas the first excited state shown in panel (b) is an 
odd function of the impurity location. The second ex- 
cited state is shown in panel (c) and this is even in the 



impurity degrees of freedom but odd in the bosonic de- 
grees of freedom specified by Nr, the number of bosons 
in the right well (the parity of the bosonic degrees of free- 
dom is defined with respect to AN). Finally, the third 
excited state is shown in panel (d) and this is odd in both 
the impurity and bosonic degrees of freedom. 



In Fig. [9] we show just the ground state, but this time 
for two different non-zero values of W; panel (a) has 
W/U = -10 and panel (b) has W/U = 10. In both cases 
the wave function has been split by the interaction with 
the impurity, but notice that the circles and the triangles 
swap positions between the two panels. In the repulsive 
case the system lowers its energy by placing more bosons 
in the opposite well to the impurity. In the limit W — > oo 
the impurity will occupy one well and all the bosons the 
other, but when the tilt is zero there is nothing to decide 
between the two wells and so the system enters a bal- 
anced macroscopic superposition of being in both wells. 
We then expect the ground state to take a Schodinger 
cat form |*°) = (|1,0) + |0, N))/\/2. In the attractive 
case the system lowers its energy by placing more bosons 
in the same well as the impurity, and again the system 
is forced into a macroscopic superposition of occupying 
both wells [41] . As W — > — oo we expect the ground state 
of the system to be |*°) = (|1, N) + |0,0))\/2- We shall 
discuss Schrodinger cat states further in Section VII 



An intriguing question concerns whether there is con- 
nection between the loop bifurcation in the mean-field 
solutions and the splitting of the many-body quantum 
state? In fact there is a direct connection, as indicated 



by Fig. 10 This figure plots W c versus Wdiy, where W c 
is the mean- field prediction, given by Eq. ( 42 1 , for the 



critical value of W at which the bifurcation in the low- 
est band occurs, and Wdi P is the value of W at which 
the ground state probability distribution for Nr first de- 
velops a dip at Nr = N/2, heralding the onset of the 
splitting. Each point in Fig. 10 corresponds to a dif- 
ferent value of the hopping energy J = J°, which lies 
in the range 1 < J/U < 180, 000. Performing a lin- 



ear fit to the data points shown in Fig. 10 we find that 
W c = —3.8 + S.lWdipi indicating that the many-body 
wave function begins splitting before the mean-field bi- 
furcation. However, the splitting of the many-body wave 
function is a smooth process (for finite N) and there is no 
particular point at which we can say the wave function 
has split. We have chosen W^ip as a simple, but arbitrary 
indicator. Another choice might be when the two peaks 
of the many-body wave function are 'resolvable', i.e. sep- 
arated by the width of the individual gaussians. With 
this indicator (not shown) we find the many-body wave 
function splits after the mean-field bifurcation. Never- 



theless, it is clear from Fig. 10 that whatever the precise 
choice of indicator there is a direct correlation between 
the occurrence of the mean-field loop and the splitting of 
the many-body wave function. 



12 



35 000 
30000 




1000 2000 3000 4000 5000 6000 7000 
W dip (units of U) 

Figure 10. (Color on line) The correlation between critical 
value W c [see Eq. ( 42 1] of the boson-impurity interaction at 
which a bifurcation occurs in the lowest band, and the equiv- 
alent value Wdip at which the many-body wave function first 
develops a dip in its center (i.e. at Nr — N/2). Each data 
point corresponds to N = 100, but to a different value of 
J = J a , lying in the range 1 < J/U < 180, 000. The solid line 
is a fit to the data and is given by W c = —3.8 + S.lVKdip. 
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Figure 11. (Color online) The coherence a, and the square of 
the number fluctuations a 2 , in the ground state as a function 
of the boson hopping energy J, with N = 100. The boson- 
impurity interaction W is set to zero, as is the tilt. The dashed 
lines indicate the boundaries between the Fock, Josephson, 
and Rabi regimes. 



VI. COHERENCE VERSUS NUMBER 
FLUCTUATIONS IN THE GROUND STATE 

Following Gati and Obcrthaler [2_2] , let us calculate the 
coherence and number fluctuations of the ground state. 
The coherence is related to the visibility of the interfer- 
ence fringes which are formed if the atoms are released 
from the double well potential, and the atomic clouds al- 
lowed to ballistically expand, until they overlap and are 
imaged in a Young's double slit type experiment. For the 
bosons the coherence is defined to be 



a = (*|B|*)/JV 



(45) 



where B = b\b 



l R 



b R b^ is the bosonic hopping operator 



In terms of the amplitudes Cm r ,n r defined in Eq. (43), 
one finds that 



Nr = N/2, and so favor a number distribution which is 
centered around Nr — N/2. The above result will only 
give the visibility of the boson fringes if the state of the 
impurity is not measured (a modified result is obtained 
if the location of the impurity is known) . 

If, rather than a Young's double slit type measurement, 
an in situ measurement of the atom number in one of the 
wells is made, then a single measurement will yield a sin- 
gle number (we assume here that the resolution of the 
measurement is at the single atom level). However, re- 
peated measurements on an identically prepared system 
will yield different atom numbers according to the prob- 
ability distribution obtained by squaring the amplitudes 
plotted in Figs. [B] and [9] The width of the atom number 
probability distribution is characterised by the standard 
deviation a, which we define via 



1 

'N 



E 

M R ,N R 



+ ^*M R ,N R + l 
N R 



x y/N R (N-N R + l) 
(Q),jv R +iCb,iVR + C* tNR _ 

x y/(N-N R )(N R + l). 



- 22 \ C Mr,n r \ 2 N r - \ ^2 \Cm r m r \ 2 N r 

M R ,N R \M R ,N R 

N R 

/ \ 2 



lCl,N R ) 



|Ci 



Nn 



Nr 



(46) 



(47) 



The coherence is therefore large when neighboring ampli- 
tudes in Fock space are large, i.e. when the number dis- 
tribution is made up of a single broad peak. On the con- 
trary, it is small if the number distribution consists of in- 
dividual spikes which are separated by at least 8Nr = 2. 
The two square root factors are peaked just either side of 



In Fig. 11 we plot the coherence and the number fluc- 
tuations for the case when the impurity is decoupled from 
the bosons (W = 0). This picture essentially repeats Fig. 
1 of [22] . We have indicated with dashed lines the bound- 
aries between three commonly defined regimes |20j : 

1. Rabi regime J/U > 2N 
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2. Josephson regime 2N > J/U > 2/N 

3. Fock regime J/U < 2/JV. 

In the Fock regime the interaction energy dominates the 
hopping energy and so the number fluctuations are small 
and the coherence low. In the Rabi regime the oppo- 
site is true, and the Josephson regime sits between the 
two. For example, Fig. |] has J/U = 10, and N = 100, 
which means that this figure corresponds to the Joseph- 
son regime. Note that when J — > oo, i.e. the d eep Rabi 
regime, then a — > \/N /2, as can be seen in Fig. [ill 

We now switch on the BEC-impurity interaction and 
ask whether a single impurity can affect the coherence 
properties of the bosons? This question is motivated by 
several previous theoretical studies [2H S3] . In [43] the 
authors consider the dynamics of a system comprising 
of two initially independent BECs which are then put 
into direct contact. They predicted the surprising result 
that macroscopic phase coherence is established between 
the two BECs as soon as a few atoms are exchanged. 
Likewise, in [24] a system comprising an atomic quantum 
dot sandwiched between two weakly coupled BECs was 
analyzed, and it was found that the dot was capable of 
mediating large amplitude Josephson-like oscillations of 
the particle imbalance despite the fact that the dot could 
only host a single atom at a time. 

The BEC-impurity system is an ideal one in which to 
study these coherence-generating mechanisms because, in 
principle, the tunnelling properties of the impurity and 
the bosons can be separately controlled. For example, 
one can examine what happens when J is very small but 
J a is not, so that the bosons are essentially frozen and 
only the impurity is mobile, thereby isolating its effect 
on the coherence between the two wells. The results are 
plotted in Fig. [l2j In this figure the value of the boson 
hopping energy is J/U = 0.001, and N — 100, which puts 
the bosons in the Fock regime, and we have set W/U = 2. 
When J a — the coherence is only a — 0.1, but as J a 
is increased it peaks at over six times this value, with 
a maximum value of a ~ 0.66 at J a /U ~ 1.5, which 
is indicated by the arrow labelled by (b) on the figure, 
before decaying back to the background value of a = 0.1 
again as J a is increased further. Note that, according to 



Eq. ( 42 ) , the mean-field prediction for the value of J a at 
which the bifurcation to form a loop in the lowest band 
occurs (when W c = 2 U) is J a = 2 U, which is close to 
the peak in a. 

To help us interpret this result we have plotted in 
Fig. [13] the many-body wave functions corresponding to 
the points labelled (a), (b), and (c) in Fig. 
virtue of the very small value of J/U = 0.001, 



t 
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le width 



(JiV/2[/) 1/4 ps 0.5 of the bosonic many-body wave func- 
tion (for either impurity state) is small and we have only 
plotted a limited range of Nr for clarity. The narrowness 
of the many-body wave function implies that the number 
fluctuations and coherence should be small, and this is 
the case except at the peak of a discussed above. Let 
us begin with panel (c), which corresponds to the case 
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Figure 12. (Color online) The coherence a, and the square of 
the number fluctuations a 2 , in the ground state as a function 
of the impurity hopping energy J a . In this figure N = 100, 
J/U=0.001, W/U = 2, and the tilt is set to zero. The three 
labels indicate the values of J a used in the different panels 
of Fig. [13] (a) J a /U = 0.001, (b) J a /U = 1.514, and (c) 
J a /U = 1000. 



where the many-body wave function is not split (and the 
mean-field loop is not present in the lowest band) be- 
cause J a /U = 1000 is so large that W c is far above the 
actual value of W = 2. The effect of the impurity is neg- 
ligible in this case, and indeed the coherence and number 
fluctuations are small. 

Jumping now to panel (a), the wave function is split 
(and the loop is present) because J a /U = 0.001 is so 
small that W c is far below the actual value of W = 2. 
The many-body wave function in this panel is to a good 
approximation given by a superposition of just two Fock 
states |*°) fa (|0,51) + |l,49))/ v / 2. Despite the fact that 
the total wave function is approximately two times wider 
than in panel (c), the coherence corresponding to (a) and 
(c) is the same, i.e. a = 0.1. Thus, even though the 
wave function in (a) has a significant total width in Fock 
space, the coherence is insensitive to this because each 
piece corresponding, respectively, to Mr = 1 and Mr = 
0, is very narrow (and we note from Eq. (46) that the 



coherence is determined separately by the two pieces of 
the wave function). Therefore, a number measurement 
is superior to a coherence measurement as a method for 
detecting the splitting because it can distinguish between 
(a) and (c). In (a) we would almost always obtain one of 
the values Nr = 49 or Nr — 51 (with equal probability), 
but almost never Nr = 50, whereas in (c) we would 
almost always obtain Nr, = 50. 

Panel (b) corresponds to the peak of a, and shows the 
wave function in the process of splitting. As can be seen, 
the parts of the wave function giving the impurity in the 
left and right wells are themselves much wider than in 
panels (a) or (c), and hence the coherence is large. In 
panel (b) we therefore have a coherent state, whereas in 
panel (c) we have (approximately) a single Fock state and 
in panel (a) we have a superposition of two Fock states. 
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Figure 13. (Color online) The many-body ground state wave 
function as a function of Nr, the number of atoms in the 
right well, for N = 100, J/U = 0.001, and W/U = 2. Each 
panel has a different value of the impurity hopping energy J a 
corresponding to the labels shown in Fig. 12 (a) J a /U = 
0.001, (b) J a /U = 1.514, and (c) J a /U = 1000. The red 
circles/blue triangles are the amplitudes for the impurity to 
be in the left/right. 



It is notable that the peak in the coherence occurs at 
the same value of J a as a step change in the number 
fluctuations a of the bosons (see lower panel of Fig. 12 ). 



As the ratio J a /U is increased and the impurity becomes 
mobile, it precipitates a sudden decrease in a at the wave 
function splitting point. 



VII. TOWARDS SCHRODINGER CAT STATES 

The value of the boson-impurity interaction which we 
have used to illustrate the coherence properties in the 
previous section, namely W = 2 U, is twice the magni- 
tude of the individual boson-boson interaction. Examin- 
ing the hamiltonian Q, and neglecting the hopping op- 
erators (this is valid in the combined Fock regime when 
J/U and J a /U are much smaller than unity), we find 
that the Fock states with AM = ±1 and AN = =f2 are 
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Figure 14. (Color online) The coherence a, and the square of 
the number fluctuations a 2 , in the ground state as a function 
of the impurity hopping energy J a . In this figure N = 100, 
J/U = 0.001, W/U = 10, and the tilt is set to zero. The 
three labels indicate the values of J a used in the different 
panels of Fig. pj (a) J a /U = 0.001, (b) J a /U = 48, and (c) 
J a /U — 1000. i\ zoom-in of this figure is given in Fig. |16| 



lower in energy than the Fock state with AM = AN = 0. 
Increasing |A7Y| beyond |A7Y| = 2 confers no further en- 
ergetic advantage when W = 2 U. This is confirmed by 
Fig. 13 (a), which is in the Fock regime, where we see that 
the two <5-function-like spikes are separated by a distance 
AN = 2. 



The state shown in Fig. 13 (a) can hardly be called a 
Schrodinger cat state because it is not a superposition 
of two macroscopically distinguishable states. Rather, it 
is a superposition of two states differing by a single bo- 
son. Indeed, this situation is rather similar to the case 
of an odd number of bosons (but no impurity) , where in 
the Fock regime the number distribution must also be 
split into two pieces differing by a single particle. In or- 
der to generate Schrodinger cat states we must increase 
the value of W: Figs. [14] and [15] illustrate the situa- 
tion when W = 10 U. We sec from Fig. 15 (a) that 
in the Fock regime we have a ground state of the form 
|*°) w (|0, 55) + |1,45))/V2, where two <5-function-like 
spikes are separated by a distance A/Vr = 10 in Fock 
space. Although perhaps not a fully fledged Schrodinger 
cat state, it is a step in that direction. 

If instead one sets W = 100 U, then one could in the- 
ory achieve the "NOON" state gU |*°) w (|0, 100) + 
|1, 0))/v^2- However, such a state would be extremely 
delicate and would quickly succumb to decoherence. In- 
terestingly, such a state is also difficult to compute be- 
cause its energy difference from the first excited state is 
exponentially small and numerical diagonalization rou- 
tines find them hard to distinguish. This situation is very 
similar to that found when computing macroscopic self- 
trapping in double- well systems (no impurity) [271 145] , 
where split number distributions are also formed. The 
problem can be diagnosed by checking if the parity of 



15 



O 

•a 

1 

O 
-o 

I 

>> 



0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
0.0 

0.4 
0.3 
0.2 
0.1 
0.0 It 



(a) 



(b) 



^ N K 



45 50 55 60 



0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
0.0 



1 


1 








(c) 









45 50 55 60 



Figure 15. (Color online) The many-body ground state wave 
function as a function of Nr, the number of atoms in the right 
well, for N = 100, J/U = 0.001, and W/U = 10. Each panel 
has a different value of the impurity hopping energy J a cor- 
responding to the labels shown in Fig. 
(b) J a /U = 48, and (c) J a /U = 1000~ 
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(a) J a /U = 0.001, 
The red circles/blue 



triangles are the amplitudes for the impurity to be in the 
left/right well. 



the eigenstates is maintained: the ground state should 
have even parity in the impurity location and the first 
excited state odd parity in the impurity location, like in 
Fig. [HJ due to the symmetry of the hamiltonian. Nu- 
merical errors can cause this parity to be lost when the 
states become nearly degenerate. One remedy for this 
situation is to form even and odd combinations of the 
(erroneously) computed ground and first excited states, 
so that parity is enforced. This method was applied for 
the smaller values of J a shown in Figs. [T4| and [To] 

The coherence and number fluctuations when W = 
10 U are shown in Fig. [Mj Four peaks can be distin- 
guished in the coherence, although the last one is quite 
wide and contains several features (a zoom-in of Fig. 14 
is given in Fig. 16 ). As before, the peaks in the coherence 



coincide with steps in the number fluctuations, and the 
three panels of Fig. [15] show the many-body wave func- 
tion in three different regimes. Panel (b) corresponds 
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Figure 16. (Color online) The coherence a, and the square of 
the number fluctuations a 2 , in the ground state as a function 
of the impurity hopping energy J a . In this figure N = 100, 
J/U = 0.001, W/U = 10, and the tilt is set to zero. This plot 
is a zoom-in of Fig. [14] 



to the maximum coherence, which is a = 0.87, and is 
located at J° = 48. The wave function at this point 
is remarkably wide considering the very small value of 
J = 0.001 U. The mean-field prediction for the bifur- 
cation point is J a = 50 U, which, once again, gives a 
good estimate for the parameter values at which dra- 
matic changes in the many-body wave function occur. 



VIII. ENTANGLEMENT ENTROPY 

One measure of the degree of entanglement between 
two subsystems is provided by their entanglement en- 
tropy. This can be calculated from the reduced density 
matrix for either subsystem. The reduced density matrix 
for the bosons is defined to be 

p B = {M R \p\M R )= (Mr\iI>)(iI>\M r ) (48) 



=0,1 



where 



|*}= E Cm r ,n r \M r ,N r ) 

M R ,N R 



(49) 



is the many-body wave function in the joint number basis 
\M R , N R ) = \M R ) \N R ). One finds that 

Pb = E (Co tN > R C 0tNR + Cl N , R Ci <Nl ^ \N R ){N' R \ . 

N R> N ' R 

.( 50 ) 

The matrix elements of the reduced density matrix for 
the bosons in the number basis are thus 



(N R \p B \N R ) - Cq n , ■ C 0j at b + C* N , ■ Ci <Nr 



(51) 



The entanglement entropy S R is defined in terms of the 
reduced density matrix as 

Sb = -J2( N r\pb^Pb\Nr) = -^p Bi lnp Bi (52) 

N R i 
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Figure 17. (Color online) Entanglement entropy S as a func- 
tion of W for N = 100 and J = J a = 1.5 U. The solid orange 
curve is for zero tilt Ae = Ae a = 0, and the dashed blue curve 
is for finite tilt Ae = Ae a = U. 



where ps i are the eigenvalues of pg. The entanglement 
entropy of the bosons is the same as that of the equivalent 
quantity for the impurity 



Si = ~y] pi t In p h 



(53) 



where the matrix elements of the reduced density matrix 
for the impurity are 



(M R \pi\M' R ) = J2 C t 



c 



Mr.Ne 



(54) 



Because of the equality of Si and Sb, we shall refer to 
the entanglement entropy simply as S rather than Sb or 
Si. 



In Fig. 17 we plot the entanglement entropy of the 
ground state as a function of W for N = 100. The solid 
orange curve is for a perfectly balanced double well and 
the dashed blue curve is for a tilted double well with 
Ae = Ae a = U. In both cases the entropy is zero 
at W — 0, because there is no entanglement when the 
boson-impurity interaction vanishes, but the two curves 
display very different behavior for large \ W\. For the bal- 
anced case the entanglement entropy approaches S = In 2 
as W —> ±oo because in that regime the bosonic ground 
state probability distribution in Fock space tends to two 
(5-function peaks at Nr — and Nr = 100 for both the 
repulsive and attractive cases. Therefore, the system can 
occupy two states each having a probability of 0.5 and 
the entropy becomes S = —2 x 0.5 In 0.5 = In 2. In the 
tilted case the ground state tends to a single (5-function 
peak in bosonic Fock space as W — > ±oo and therefore 
the entanglement entropy tends to zero in those limits. 
Of course, the tilt can never be precisely zero in an ex- 
periment and so in reality we always expect S — > if the 
limit W — > ±oo can be achieved. However, as the tilt be- 
comes small the value of \W\ required for the system to 
become sensitive to a finite value of Ae or Ae a becomes 
large. 



Once again, it is interesting to ask whether the mean- 
field prediction for the bifurcation point given by Eq. ( 42 ) 
bears any relevance to the many-body results. For the 
parameters used in Fig. [17J we find that W c = ±1.76 U 
(this result is only for the zero tilt case). These values 
do not correspond to any obvious features on the solid 
orange curve in Fig. |17| but we note that the first deriva- 
tive of this curve with respect to W has extrema close by 
at W = ±2.68 l7. 



IX. SUMMARY AND DISCUSSION 

In this paper we have studied the effects a single impu- 
rity can have upon a BEC in a double well potential, with 
the emphasis placed upon reconciling the many-body and 
mean-field descriptions. The mean-field theory is nonlin- 
ear (indeed, our system is analogous to a double pendu- 
lum, which is known to be chaotic) and the static solu- 
tions undergo bifurcations as parameters such as the tilt 
Ae, and the boson-impurity interaction W, are varied. 
In particular, the lowest lying mean-field solution (cor- 
responding to both pendula pointing down) undergoes a 
pitchfork bifurcation when W is increased past a critical 



value W c given by Eq. ( 42 ) . This bifurcation is due to the 



spontaneous formation of an imbalance in the number of 
atoms in the two wells. Considered as a function of the 
tilt between the two wells, the bifurcated solutions form 
a swallowtail loop in the lowest lying mean-field solution. 
This critical value of the boson-impurity interaction need 
not be large: taking numbers from the experiments [3] 
and [5] and extrapolating them to include an impurity 
allows one to estimate that W c ~ 10 U. 

The bifurcation in the lowest lying mean-field solution 
corresponds, in the many-body theory, to a splitting of 
the ground state atom number probability distribution 
into two separate pieces. One piece corresponds to the 
impurity being in the left well and an increased (for posi- 
tive values of W) or decreased (for negative values of W) 
number of bosons in the right well, and vice versa for the 
other piece. Unlike the mean-field solution, the many- 
body solution is in general a superposition of both cases 
(at least until decoherence is added into the model), and 
in the large W limit takes the form of a Schrodinger cat 
state. 

The presence of the impurity can have a dramatic ef- 
fect on the coherence a of the bosons between the two 
wells: at the bifurcation a is strongly peaked, even for 
exceedingly small values of the bosonic hopping energy 
J. This phenomenon provides a readily verifiable ex- 
perimental signature of the presence of the bifurcation. 
A single particle tunneling between the two wells can 
therefore drive the system into a coherent state [52] 33] . 
However, either side of the bifurcation a returns to the 
very low background level set by J, and, surprisingly, 
this is true even when J a is very large so that the im- 
purity is highly mobile. It is conceivable that the high 
coherence at the bifurcation in a BEC-impurity system 
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could be usefully employed in interferometry [T2], where 
one wants both high coherence and yet also small number 
fluctuations (because the boson-boson interactions mean 
that different boson numbers in the wells lead to different 
mean-field shifts). 

From the quantum measurement theory perspective, 
the formation of a Schrodinger cat state is necessary if the 
BEC is to act as a macroscopic device that measures the 
position of the microscopic impurity. For, only then can 
the BEC unambiguously distinguish between the impu- 
rity being in the left and right wells (the Schrodinger cat 
state is assumed to collapse onto one of its two compos- 
ite states under the influence of an environment). Given 
that the bifurcation heralds the splitting of the number 
distribution, it is interesting to inquire about the fate of 
the bifurcation in the macroscopic (large N) limit. Ex- 
amining Eq. (42 1, we see that if NU J , then W c tends 
to the value 



lim WAN) 



V2JHJ 



(55) 



which is independent of N. This expression bears a re- 
semblance to Eq. ( |33[ ) for the plasma frequency, which 
gives the energy of the first excited state of a BEC in 
a double well in the same limit, except that J has been 
replaced by J a and N has been set to unity. 



Taken at face value, Eq. (55) seems to imply the non- 



sensical result that a single impurity with a finite inter- 
action W can have a finite effect on an infinitely large 
system. However, this interpretation is misleading. In 
particular, Eq. (55) is written in terms of microscopic 



parameters and should instead be expressed in terms of 
intensive quantities that have meaning in the thermody- 
namic limit N —J- oo, V —> oo, but N/V ^constant. To 
accomplish this we note that W is the interaction energy 
per particle, and should properly be compared with the 
plasma energy per particle /iWpi as /iV. This gives 



lim WAN) 

N-y+co 



Hui 



plas 



N 



N . 



(56) 



where we have put J = J a . Taking the plasma energy 
per particle as an intensive quantity that is independent 
of N, we find that W c scales as y/~N as N — > oo, and 
thus an infinite boson-impurity interaction is required to 
trigger a bifurcation in the thermodynamic limit. 
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Appendix A: The double pendulum 

For the convenience of the reader, in this Appendix 
we summarize some results concerning the double pen- 
dulum, which is a system made up of one pendulum sus- 
pended from another. Each pendulum consists of a mass- 
less rod of length ij, a bob of mass m^, and subtends an 
angle 9i to the downward vertical. The upper pendulum 
corresponds to i = 1 and the lower pendulum to i = 2. 
The kinetic T and potential V energies for this system 
are 0S1 



T 



m,\l\9\ m 2 
2 + IT 



l\9\ + l\Q\ + 2l x l 2 9 x 9 2 cos(9 1 - 9 2 ) 



V = —(mi + m 2 )gl\ cos6q — m 2 gl 2 cos #2 



(Al) 
(A2) 



In the main part of the text, we describe the BEC- 
impurity system in terms of the phase angles a and /3, 
and their conjugate number differences Y and Z. In or- 
der to express the double pendulum in terms of conjugate 
variables, we form the lagrangian L = T — V , and obtain 
the conjugate momenta via Pi = dL/dOi. We find 



pi = (mi + m 2 )l\9\ + m 2 lil 2 9 2 cos(#i - 9 2 ) 



p 2 = m 2 l\9 2 



m 2 l 1 l 2 9 1 cos(6»! - 2 ) 



(A3) 
(A4) 



Solving for 9i and 9 2l we can eliminate these angular 
velocities from the hamiltonian H = Bipi — L in favor of 
the conjugate momenta to give 



H, 



dp 



llp 2 (m 1 + m 2 ) + l\p\m 2 - 2lil 2 pip 2 m 2 cos(#i - 9 2 



2l\l 2 m 2 \m\ + m 2 sin (6>i - 
-(mi + m 2 )gl\ cos Q\ — m 2 gl 2 cos 9 2 



(A5) 



We have some freedom to choose the two masses since, 
when uncoupled, the frequencies of the two pendula do 
not depend upon them. In particular, when mi 3> m,2 
we have 



H, 



dp 



1 r pi 

2 \ mil\ 



_P|_ 
m 2 l\ 



2 Pi Pi 
2 mi h h 
m\gl\ cos 9\ — m 2 gl 2 cos 9 2 . 
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(A6) 



For example, we might take the more massive pendulum 
to correspond to the BEC and the less massive one to the 
impurity. Providing the two pendula are close to station- 
ary points, i.e. where both of them are either pointing 
downwards or upwards or one is pointing upwards and 
the other downwards, the cos(#i — 9 2 ) term can be re- 
placed with ±1, as appropriate. This allows for a close 
correspondence with the mean-field hamiltonian Q8J) , al- 
though the sign of the boson-impurity interaction W is 
then set depending upon whether cos(#i — 9 2 ) equals +1 
or —1. 

Although it might appear that the term corresponding 
to p 2 , i.e. a term in Y 2 , is missing from the hamiltonian 
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(251, this is not the case providing Y <C 1. For then we can expand the \/T — 4F 2 term to give 

H w C/Z 2 - J^/N 2 - 4Z 2 cos /3 + 2 J Q F 2 - J a cos a 
+2WYZ . (A7) 

where we have put Ae = Ae a = 0. Furthermore, putting 
y/N 2 - 4Z 2 -> N, which assumes NU ^> J, we obtain a 
hamiltonian equivalent to Eq. ( A6 1 . 
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